using System;
using System.Collections.Generic;
using System.Text;

namespace ImagesMath
{
    public class M33
    {

        private VR3 vx, vy, vz;

        #region Constructors

        public M33()
        {
            double a =1;
            double b =0;
            double c =0;
            double d =0;
            double e =1;
            double f =0;
            double g =0;
            double h =0;
            double i =1;
            bool column = true;

            if (column)
            {

                vx = new VR3(a, b, c);
                vy = new VR3(d, e, f);
                vz = new VR3(g, h, i);
            }
            else
            {
                vx = new VR3(a, d, g);
                vy = new VR3(b, e, h);
                vz = new VR3(c, f, i);
            }
            //if (column)
            //{
            //    VX.X = a; VY.X = d; VZ.X = g;
            //    VX.Y = b; VY.Y = e; VZ.Y = h;
            //    VX.Z = c; VY.Z = f; VZ.Z = i;
            //}
            //else
            //{
            //    VX.X = a; VY.X = b; VZ.X = c;
            //    VX.Y = d; VY.Y = e; VZ.Y = f;
            //    VX.Z = g; VY.Z = h; VZ.Z = i;
            //}
        }

        public M33( double a,double b,double c,
                    double d,double e,double f,
                    double g,double h,double i,bool column)
        {
            if (column)
            {

                vx = new VR3(a, b, c);
                vy = new VR3(d, e, f);
                vz = new VR3(g, h, i);
            }
            else
            {
                vx = new VR3(a, d, g);
                vy = new VR3(b, e, h);
                vz = new VR3(c, f, i);
            }
        }
        public M33(VR3 x,VR3 y,VR3 z,bool column)
        {
            if (column)
            {
                vx = x; vy = y; vz = z;
            }
            else
            {
                vx = new VR3(x.X, x.Y, x.Z);
                vy = new VR3(y.X, y.Y, y.Z);
                vz = new VR3(z.X, z.Y, z.Z);
                
                //VX.X = x.X; VX.Y = y.X; VX.Z = z.X;
                //VY.X = x.Y; VY.Y = y.Y; VY.Z = z.Y;
                //VZ.X = x.Z; VZ.Y = y.Z; VZ.Z = z.Z;
            }
        }

        #endregion

        #region Accessors

        public VR3 VX
        {
            get { return this.vx; }
            set { this.vx = value; }
        }
        public VR3 VY
        {
            get { return this.vy; }
            set { this.vy = value; }
        }
        public VR3 VZ
        {
            get { return this.vz; }
            set { this.vz = value; }
        }

        #endregion

        #region Arithmetical operators

        public static M33 operator + (M33 m1, M33 m2 )
        {
            return new M33(m1.VX + m2.VX, m1.VY + m2.VY, m1.VZ + m2.VZ, true);
        }


        public static M33 operator -(M33 m1, M33 m2)
        {
            return new M33(m1.VX - m2.VX, m1.VY - m2.VY, m1.VZ - m2.VZ, true);
        }


        public static M33 operator -(M33 m1, double k)
        {
            return new M33(m1.VX * k, m1.VY * k, m1.VZ * k, true);
        }


        //m33 & m33::operator *= (double k)
        //{
        //  vx*=k; vy*=k; vz*=k;
        //  return *this;
        //}

        public static PR3 operator *(M33 m, PR3 p)
        {
            return new PR3( p.X * m.VX.X + p.Y * m.VY.X + p.Z * m.VZ.X,
                            p.X * m.VX.Y + p.Y * m.VY.Y + p.Z * m.VZ.Y,
                            p.X * m.VX.Z + p.Y * m.VY.Z + p.Z * m.VZ.Z);
        }

        //pR3 m33::operator * (const pR3 &p) const
        //{
        //  return pR3(p.x*VX.X + p.y*VY.X + p.z*VZ.X,
        //    p.x*VX.Y + p.y*VY.Y + p.z*VZ.Y,
        //    p.x*VX.Z + p.y*VY.Z + p.z*VZ.Z);
        //}

        public static VR3 operator *(M33 m, VR3 v)
        {
            return new VR3( v.X * (m.VX).X + v.Y * m.VY.X + v.Z * m.VZ.X,
                            v.X * m.VX.Y + v.Y * m.VY.Y + v.Z * m.VZ.Y,
                            v.X * m.VX.Z + v.Y * m.VY.Z + v.Z * m.VZ.Z);
        }

        //vR3 m33::operator * (const vR3 &u) const
        //{
        //  return vR3(u.x*VX.X + u.y*VY.X + u.z*VZ.X,
        //    u.x*VX.Y + u.y*VY.Y + u.z*VZ.Y,
        //    u.x*VX.Z + u.y*VY.Z + u.z*VZ.Z);
        //}


        public static M33 operator *(M33 m1, M33 m2)
        {
            return new M33( m1.VX.X * m2.VX.X + m1.VY.X * m2.VX.Y + m1.VZ.X * m2.VX.Z,
                            m1.VX.Y * m2.VX.X + m1.VY.Y * m2.VX.Y + m1.VZ.Y * m2.VX.Z,
                            m1.VX.Z * m2.VX.X + m1.VY.Z * m2.VX.Y + m1.VZ.Z * m2.VX.Z,
                            m1.VX.X * m2.VY.X + m1.VY.X * m2.VY.Y + m1.VZ.X * m2.VY.Z,
                            m1.VX.Y * m2.VY.X + m1.VY.Y * m2.VY.Y + m1.VZ.Y * m2.VY.Z,
                            m1.VX.Z * m2.VY.X + m1.VY.Z * m2.VY.Y + m1.VZ.Z * m2.VY.Z,
                            m1.VX.X * m2.VZ.X + m1.VY.X * m2.VZ.Y + m1.VZ.X * m2.VZ.Z,
                            m1.VX.Y * m2.VZ.X + m1.VY.Y * m2.VZ.Y + m1.VZ.Y * m2.VZ.Z,
                            m1.VX.Z * m2.VZ.X + m1.VY.Z * m2.VZ.Y + m1.VZ.Z * m2.VZ.Z,true);
        }

        //m33 m33::operator * (const m33 &m) const
        //{
        //  return m33(VX.X*m.VX.X + VY.X*m.VX.Y + VZ.X*m.VX.Z,
        //    VX.Y*m.VX.X + VY.Y*m.VX.Y + VZ.Y*m.VX.Z,
        //    VX.Z*m.VX.X + VY.Z*m.VX.Y + VZ.Z*m.VX.Z,
        //    VX.X*m.VY.X + VY.X*m.VY.Y + VZ.X*m.VY.Z,
        //    VX.Y*m.VY.X + VY.Y*m.VY.Y + VZ.Y*m.VY.Z,
        //    VX.Z*m.VY.X + VY.Z*m.VY.Y + VZ.Z*m.VY.Z,
        //    VX.X*m.VZ.X + VY.X*m.VZ.Y + VZ.X*m.VZ.Z,
        //    VX.Y*m.VZ.X + VY.Y*m.VZ.Y + VZ.Y*m.VZ.Z,
        //    VX.Z*m.VZ.X + VY.Z*m.VZ.Y + VZ.Z*m.VZ.Z);
        //}


        //m33 & m33::operator *= (const m33 &m)
        //{
        //  *this=*this*m;
        //  return *this;
        //}

        public static double operator ! (M33 m ) 
        {
            return (m.VX.X * m.VX.X + m.VX.Y * m.VX.Y + m.VX.Z * m.VX.Z +
            m.VY.X * m.VY.X + m.VY.Y * m.VY.Y + m.VY.Z * m.VY.Z +
            m.VZ.X * m.VZ.X + m.VZ.Y * m.VZ.Y + m.VZ.Z * m.VZ.Z);
        }

        //double m33::operator ! ( void ) const
        //{
        //  return (VX.X*VX.X + VX.Y*VX.Y + VX.Z*VX.Z +
        //    VY.X*VY.X + VY.Y*VY.Y + VY.Z*VY.Z +
        //    VZ.X*VZ.X + VZ.Y*VZ.Y + VZ.Z*VZ.Z);
        //}

        #endregion

        #region Functions

        public M33 transpose()
        {
          return new M33(vx,vy,vz,false);
        }


        //void m33::tridiagonalize(double Diag[3], double SubD[3])
        //{
        ///// Householder reduction T = Q^t M Q
        /////   Input:
        ///// 	 mat, symmetric 3x3 matrix M
        /////   Output:
        ///// 	 mat, orthogonal matrix Q (a reflection)
        ///// 	 diag, diagonal entries of T
        ///// 	 subd, subdiagonal entries of T (T is symmetric)

        //  double fA(VX.X), fB(VY.X), fC(VZ.X), fD(VY.Y), fE(VZ.Y), fF(VZ.Z);

        //  Diag[0] = fA;
        //  SubD[2] = 0.0;
        //  if (fabs(fC) >= DBL_EPSILON)
        //  {
        //    double length = sqrt(SQ(fB)+SQ(fC));
        //    double invlen = 1.0/length;
        //    fB *= invlen;
        //    fC *= invlen;
        //    double fQ = 2.0*fB*fE+fC*(fF-fD);
        //    Diag[1] = fD+fC*fQ;
        //    Diag[2] = fF-fC*fQ;
        //    SubD[0] = length;
        //    SubD[1] = fE-fB*fQ;

        //    vx = vR3::UNIT_X;
        //    vy = vR3(0.0,fB,fC);
        //    vz = vR3(0.0,fC,-fB);
        //  }
        //  else
        //  {
        //    Diag[1] = fD;
        //    Diag[2] = fF;
        //    SubD[0] = fB;
        //    SubD[1] = fE;

        //    vx =   vR3::UNIT_X;
        //    vy =   vR3::UNIT_Y;
        //    vz = - vR3::UNIT_Z;                           // A verif
        //  }
        //}


        //bool m33::QLAlgorithm(double Diag[3], double SubD[3])
        //{
        ///// QL iteration with implicit shifting to reduce matrix from tridiagonal
        ///// to diagonal
        //  double mtab[9],fTmp0,fTmp1,fTmp2,fTmp3,fTmp4,fCos,fSin;

        //  mtab[0] = VX.X; mtab[1] = VY.X; mtab[2] = VZ.X;
        //  mtab[3] = VX.Y; mtab[4] = VY.Y; mtab[5] = VZ.Y;
        //  mtab[6] = VX.Z; mtab[7] = VY.Z; mtab[8] = VZ.Z;

        //  for (int i0 = 0; i0 < 3; i0++)
        //  {
        //    const int iMaxIter = 32;
        //    int iIter;
        //    for (iIter = 0; iIter < iMaxIter; iIter++)
        //    {
        //      int i1;
        //      for (i1 = i0; i1 <= 1; i1++)
        //      {
        //        double fSum = fabs(Diag[i1]) + fabs(Diag[i1+1]);
        //        if (fabs(SubD[i1]) + fSum == fSum) break;
        //      }
        //      if (i1 == i0) break;

        //      fTmp0 = (Diag[i0+1] - Diag[i0]) / (2.0*SubD[i0]);
        //      fTmp1 = sqrt(SQ(fTmp0)+1.0);

        //      fTmp0 = Diag[i1]-Diag[i0]+SubD[i0]/(fTmp0+SIGN(fTmp1,fTmp0));

        //      fSin = fCos = 1.0;
        //      fTmp2 = 0.0;
        //      for (int i2 = i1-1; i2 >= i0; i2--)
        //      {
        //        fTmp3 = fSin*SubD[i2];
        //        fTmp4 = fCos*SubD[i2];

        //        if (fabs(fTmp3) >= fabs(fTmp0))
        //        {
        //          fCos = fTmp0/fTmp3;
        //          fTmp1 = sqrt(SQ(fCos)+1.0);
        //          SubD[i2+1] = fTmp3*fTmp1;
        //          fCos *= (fSin = 1.0/fTmp1);
        //        }
        //        else
        //        {
        //          fSin = fTmp3/fTmp0;
        //          fTmp1 = sqrt(SQ(fSin)+1.0);
        //          SubD[i2+1] = fTmp0*fTmp1;
        //          fSin *= (fCos = 1.0/fTmp1);
        //        }
        //        fTmp0 = Diag[i2+1]-fTmp2;
        //        fTmp1 = (Diag[i2]-fTmp0)*fSin+2.0*fCos*fTmp4;
        //        fTmp2 = fSin*fTmp1;
        //        Diag[i2+1] = fTmp0+fTmp2;
        //        fTmp0 = fCos*fTmp1-fTmp4;

        //        for (int iRow = 0; iRow < 3; iRow++)
        //        {
        //          fTmp3 = mtab[(3*iRow)+i2+1];
        //          mtab[(3*iRow)+i2+1] = fSin*mtab[(3*iRow)+i2] + fCos*fTmp3;
        //          mtab[(3*iRow)+i2]   = fCos*mtab[(3*iRow)+i2] - fSin*fTmp3;
        //        }
        //      }
        //      Diag[i0] -= fTmp2;
        //      SubD[i0]  = fTmp0;
        //      SubD[i1]  = 0.0;
        //    }
        //    assert (iIter != iMaxIter) ;
        //  }
        //  VX.X = mtab[0]; VY.X = mtab[1]; VZ.X = mtab[2];
        //  VX.Y = mtab[3]; VY.Y = mtab[4]; VZ.Y = mtab[5];
        //  VX.Z = mtab[6]; VY.Z = mtab[7]; VZ.Z = mtab[8];

        //  return true;
        //}


        //double m33::norm2( void ) const
        //{
        //  return (VX.X*VX.X + VX.Y*VX.Y + VX.Z*VX.Z +
        //    VY.X*VY.X + VY.Y*VY.Y + VY.Z*VY.Z +
        //    VZ.X*VZ.X + VZ.Y*VZ.Y + VZ.Z*VZ.Z);
        //}


        //double m33::norm( void ) const
        //{
        //  return sqrt(norm2());
        //}


        //double m33::trace() const
        //{
        //  return VX.X+VY.Y+VZ.Z;
        //}

        public double determinant()
        {
          return   VX.X * VY.Y * VZ.Z
            + VX.Y * VY.Z * VZ.X
            + VX.Z * VY.X * VZ.Y
            - VX.Z * VY.Y * VZ.X
            - VX.Y * VY.X * VZ.Z
            - VX.X * VY.Z * VZ.Y;
        }


        public M33 inverse() 
        {
          double det, pos, neg, temp;
          

        // Calculate the determinant of submatrix A
          pos = neg = 0.0f;
          temp =  VX.X * VY.Y * VZ.Z; if (temp >= 0.0f) pos += temp;  else  neg += temp;
          temp =  VX.Y * VY.Z * VZ.X; if (temp >= 0.0f) pos += temp;  else  neg += temp;
          temp =  VX.Z * VY.X * VZ.Y; if (temp >= 0.0f) pos += temp;  else  neg += temp;
          temp = -VX.Z * VY.Y * VZ.X; if (temp >= 0.0f) pos += temp;  else  neg += temp;
          temp = -VX.Y * VY.X * VZ.Z; if (temp >= 0.0f) pos += temp;  else  neg += temp;
          temp = -VX.X * VY.Z * VZ.Y; if (temp >= 0.0f) pos += temp;  else  neg += temp;
          det = pos + neg;


            

        // Is Singular ?
          //assert(!((det == 0.0f) || (fabs(det / (pos - neg)) < 1e-15)));

        // Calculate inverse(A) = adj(A) / det(A)
          det = 1.0f / det;
            VR3 vx = new VR3(( (VY.Y * VZ.Z - VY.Z * VZ.Y) * det),(-(VX.Y * VZ.Z - VX.Z * VZ.Y) * det),( (VX.Y * VY.Z - VX.Z * VY.Y) * det));

            //inv.VX.X =( (VY.Y * VZ.Z - VY.Z * VZ.Y) * det);
            //inv.VX.Y =(-(VX.Y * VZ.Z - VX.Z * VZ.Y) * det);
            //inv.VX.Z =( (VX.Y * VY.Z - VX.Z * VY.Y) * det);
          
          

            VR3 vy = new VR3((-(VY.X * VZ.Z - VY.Z * VZ.X) * det), ((VX.X * VZ.Z - VX.Z * VZ.X) * det), (-(VX.X * VY.Z - VX.Z * VY.X) * det));


            //inv.VY.X =(-(VY.X * VZ.Z - VY.Z * VZ.X) * det);
            //inv.VY.Y =( (VX.X * VZ.Z - VX.Z * VZ.X) * det);
            //inv.VY.Z =(-(VX.X * VY.Z - VX.Z * VY.X) * det);
            
          

            VR3 vz = new VR3(((VY.X * VZ.Y - VY.Y * VZ.X) * det), (-(VX.X * VZ.Y - VX.Y * VZ.X) * det), ((VX.X * VY.Y - VX.Y * VY.X) * det));

            //inv.VZ.X =( (VY.X * VZ.Y - VY.Y * VZ.X) * det);
            //inv.VZ.Y =(-(VX.X * VZ.Y - VX.Y * VZ.X) * det);
            //inv.VZ.Z =( (VX.X * VY.Y - VX.Y * VY.X) * det);

            M33 inv = new M33(vx,vy,vz,true);

          return inv;
        }

        //m33 m33::inverse() const
        //{
        //  double det, pos, neg, temp;
        //  m33 inv;

        //// Calculate the determinant of submatrix A
        //  pos = neg = 0.0f;
        //  temp =  VX.X * VY.Y * VZ.Z; if (temp >= 0.0f) pos += temp;  else  neg += temp;
        //  temp =  VX.Y * VY.Z * VZ.X; if (temp >= 0.0f) pos += temp;  else  neg += temp;
        //  temp =  VX.Z * VY.X * VZ.Y; if (temp >= 0.0f) pos += temp;  else  neg += temp;
        //  temp = -VX.Z * VY.Y * VZ.X; if (temp >= 0.0f) pos += temp;  else  neg += temp;
        //  temp = -VX.Y * VY.X * VZ.Z; if (temp >= 0.0f) pos += temp;  else  neg += temp;
        //  temp = -VX.X * VY.Z * VZ.Y; if (temp >= 0.0f) pos += temp;  else  neg += temp;
        //  det = pos + neg;

        //// Is Singular ?
        //  assert(!((det == 0.0f) || (fabs(det / (pos - neg)) < 1e-15)));

        //// Calculate inverse(A) = adj(A) / det(A)
        //  det = 1.0f / det;
        //  inv.VX.X =( (VY.Y * VZ.Z - VY.Z * VZ.Y) * det);
        //  inv.VY.X =(-(VY.X * VZ.Z - VY.Z * VZ.X) * det);
        //  inv.VZ.X =( (VY.X * VZ.Y - VY.Y * VZ.X) * det);
        //  inv.VX.Y =(-(VX.Y * VZ.Z - VX.Z * VZ.Y) * det);
        //  inv.VY.Y =( (VX.X * VZ.Z - VX.Z * VZ.X) * det);
        //  inv.VZ.Y =(-(VX.X * VZ.Y - VX.Y * VZ.X) * det);
        //  inv.VX.Z =( (VX.Y * VY.Z - VX.Z * VY.Y) * det);
        //  inv.VY.Z =(-(VX.X * VY.Z - VX.Z * VY.X) * det);
        //  inv.VZ.Z =( (VX.X * VY.Y - VX.Y * VY.X) * det);

        //  return inv;
        //}


        //m33 m33::transpose() const
        //{
        //  return m33(vx,vy,vz,false);
        //}


        //m33& m33::orthonormalize()
        //{
        ///// Algorithm uses Gram-Schmidt orthogonalization.  If 'this' matrix is
        ///// M = [m0|m1|m2], then orthonormal output matrix is Q = [q0|q1|q2],
        /////
        /////   q0 = m0/|m0|
        /////   q1 = (m1-(q0*m1)q0)/|m1-(q0*m1)q0|
        /////   q2 = (m2-(q0*m2)q0-(q1*m2)q1)/|m2-(q0*m2)q0-(q1*m2)q1|
        /////
        ///// where |V| indicates length of vector V and A*B indicates dot
        ///// product of vectors A and B.

        //// compute q0
        //  vx.normalize();

        //// compute q1
        //  double dot0(vx%vy);
        //  vy -= vx*dot0;
        //  vy.normalize();

        //// compute q2
        //  double dot1(vy%vz);
        //  dot0 = vx%vz;

        //  vz -= vx*dot0 + vy*dot1;
        //  vz.normalize();

        //  return *this;
        //}


        //void m33::fromEulerAnglesXYZ(double YAngle, double PAngle, double RAngle)
        //{
        //  double c1(cos(YAngle)), s1(sin(YAngle));
        //  double c2(cos(PAngle)), s2(sin(PAngle));
        //  double c3(cos(RAngle)), s3(sin(RAngle));
        //  m33 xmat(1.0,0.0,0.0,0.0,c1,-s1,0.0,s1,c1,false);
        //  m33 ymat(c2,0.0,s2,0.0,1.0,0.0,-s2,0.0,c2,false);
        //  m33 zmat(c3,-s3,0.0,s3,c3,0.0,0.0,0.0,1.0,false);

        //  *this = xmat*(ymat*zmat);
        //}


        //bool m33::toEulerAnglesXYZ(double &XAngle, double &YAngle, double &ZAngle) const
        //{
        //// rot =  cy*cz			 -cy*sz  			sy
        //// 		 cz*sx*sy+cx*sz  cx*cz-sx*sy*sz -cy*sx
        //// 		-cx*cz*sy+sx*sz  cz*sx+cx*sy*sz  cx*cy

        //  if ( VZ.X < 1.0 )
        //  {
        //    if ( VZ.X > -1.0 )
        //    {
        //      XAngle = atan2(-VZ.Y,VZ.Z);
        //      YAngle = asin((double)VZ.X);
        //      ZAngle = atan2(-VY.X,VX.X);
        //      return true;
        //    }
        //    else
        //    {
        //// WARNING. Not unique
        //      XAngle = -atan2(VX.Y,VY.Y);
        //      YAngle = -HALF_PI;
        //      ZAngle = 0.0;
        //      return false;
        //    }
        //  }
        //  else
        //  {
        //// WARNING. Not unique
        //    XAngle = atan2(VX.Y,VY.Y);
        //    YAngle = HALF_PI;
        //    ZAngle = 0.0;
        //    return false;
        //  }
        //}


        //void m33::toAxisAngle(vR3 &axe, double &angle) const
        //{
        ///// Let (x,y,z) be the unit-length axis and let A be an angle of rotation.
        ///// The rotation matrix is R = I + sin(A)*P + (1-cos(A))*P^2 where
        ///// I is the identity and
        /////
        ///// 		+- 		 -+
        ///// P = |  0 -z +y |
        ///// 		| +z  0 -x |
        ///// 		| -y +x  0 |
        ///// 		+- 		 -+
        /////
        ///// If A > 0, R represents a counterclockwise rotation about the axis in
        ///// the sense of looking from the tip of the axis vector towards the
        ///// origin.  Some algebra will show that
        /////
        /////   cos(A) = (trace(R)-1)/2  and  R - R^t = 2*sin(A)*P
        /////
        ///// In the event that A = pi, R-R^t = 0 which prevents us from extracting
        ///// the axis through P.  Instead note that R = I+2*P^2 when A = pi, so
        ///// P^2 = (R-I)/2.  The diagonal entries of P^2 are x^2-1, y^2-1, and
        ///// z^2-1.  We can solve these for axis (x,y,z).  Because the angle is pi,
        ///// it does not matter which sign you choose on the square roots.

        //  double tr(trace()), cs(0.5*(tr-1.0));

        //  angle = acos(cs);

        //  if (angle > 0.0)
        //  {
        //    if (angle < PI)
        //    {
        //      axe.x = VY.Z-VZ.Y;
        //      axe.y = VZ.X-VX.Z;
        //      axe.z = VX.Y-VY.X;
        //      axe.normalize();
        //    }
        //    else
        //    {
        //      double hInv;
        //      if (VX.X >= VY.Y)
        //      {
        //// m00 >= m11
        //        if (VX.X >= VZ.Z)
        //        {
        //// m00 is maximum diagonal term
        //          axe.x = 0.5*sqrt(VX.X - VY.Y - VZ.Z + 1.0);
        //          hInv = 0.5/axe.x;
        //          axe.y = hInv*VY.X;
        //          axe.z = hInv*VZ.X;
        //        }
        //        else
        //        {
        //// m22 is maximum diagonal term
        //          axe.z = 0.5*sqrt(VZ.Z - VX.X - VY.Y + 1.0);
        //          hInv = 0.5/axe.z;
        //          axe.x = hInv*VZ.X;
        //          axe.y = hInv*VZ.Y;
        //        }
        //      }
        //      else
        //      {
        //// m11 > m00
        //        if (VY.Y >= VZ.Z)
        //        {
        //// m11 is maximum diagonal term
        //          axe.y = 0.5*sqrt(VY.Y - VX.X - VZ.Z + 1.0);
        //          hInv  = 0.5/axe.y;
        //          axe.x = hInv*VY.X;
        //          axe.z = hInv*VZ.Y;
        //        }
        //        else
        //        {
        //// m22 is maximum diagonal term
        //          axe.z = 0.5*sqrt(VZ.Z - VX.X - VY.Y + 1.0);
        //          hInv = 0.5/axe.z;
        //          axe.x = hInv*VZ.X;
        //          axe.y = hInv*VZ.Y;
        //        }
        //      }
        //    }
        //  }
        //// The angle is 0 and the matrix is the identity.  Any axis will
        //// work, so just use the x-axis.
        //  else axe = vR3::UNIT_X;
        //}


        //void m33::fromAxisAngle(const vR3 &axe, double angle)
        //{
        //  double fCos(cos(angle)), fSin(sin(angle)), fOneMinusCos(1.0-fCos);
        //  double fX2(SQ(axe.x)), fY2(SQ(axe.y)), fZ2(SQ(axe.z));
        //  double fXYM(axe.x*axe.y*fOneMinusCos);
        //  double fXZM(axe.x*axe.z*fOneMinusCos);
        //  double fYZM(axe.y*axe.z*fOneMinusCos);
        //  double fXSin(axe.x*fSin), fYSin(axe.y*fSin), fZSin(axe.z*fSin);

        //  VX.X = fX2*fOneMinusCos+fCos;
        //  VX.Y = fXYM+fZSin;
        //  VX.Z = fXZM-fYSin;
        //  VY.X = fXYM-fZSin;
        //  VY.Y = fY2*fOneMinusCos+fCos;
        //  VY.Z = fYZM+fXSin;
        //  VZ.X = fXZM+fYSin;
        //  VZ.Y = fYZM-fXSin;
        //  VZ.Z = fZ2*fOneMinusCos+fCos;
        //}


        //void m33::toGL(double *m) const
        //{
        //  m[0]  = VX.X;
        //  m[1]  = VX.Y;
        //  m[2]  = VX.Z;
        //  m[3]  = 0.;

        //  m[4]  = VY.X;
        //  m[5]  = VY.Y;
        //  m[6]  = VY.Z;
        //  m[7]  = 0.;

        //  m[8]  = VZ.X;
        //  m[9]  = VZ.Y;
        //  m[10] = VZ.Z;
        //  m[11] = 0.;

        //  m[12] = 0;
        //  m[13] = 0;
        //  m[14] = 0;
        //  m[15] = 1.;
        //}


        //m33 m33::slerp(double t, const m33 &m0, const m33 &m1)
        //{
        //  double angle;
        //  vR3 axe;
        //  m33 m(m0.transpose()*m1);

        //  m.toAxisAngle(axe,angle);
        //  m.fromAxisAngle(axe,t*angle);

        //  return m;
        //}


        //void m33::eigenDecomposition(m33 &R, m33 &D) const
        //{
        //// Factor M = R*D*R^T.  The columns of R are the eigenvectors.  The
        //// diagonal entries of D are the corresponding eigenvalues.
        //  double save, diag[3], subdiag[3];

        //  R = *this;
        //  R.tridiagonalize(diag,subdiag);
        //  R.QLAlgorithm   (diag,subdiag);

        //// The Householder transformation is a reflection.  Make the eigenvectors
        //// a right--handed system by changing sign on the last column.

        ////R.VX.Z = -R.VX.Z; R.VY.Z = -R.VY.Z; R.VZ.Z = -R.VZ.Z; // A verif

        //// sort eigenvalues in increasing order, d0 <= d1 <= d2
        //  if (diag[1] < diag[0])
        //  {
        //// swap d0 and d1
        //    save = diag[0]; diag[0] = diag[1]; diag[1] = save;

        //// swap V0 and V1
        //    save = R.VX.X; R.VX.X = R.VY.X; R.VY.X = save;
        //    save = R.VX.Y; R.VX.Y = R.VY.Y; R.VY.Y = save;
        //    save = R.VX.Z; R.VX.Z = R.VY.Z; R.VY.Z = save;
        //  }

        //  if (diag[2] < diag[1])
        //  {
        //// swap d1 and d2
        //    save = diag[1]; diag[1] = diag[2]; diag[2] = save;

        //// swap V1 and V2
        //    save = R.VY.X; R.VY.X = R.VZ.X; R.VZ.X = save;
        //    save = R.VY.Y; R.VY.Y = R.VZ.Y; R.VZ.Y = save;
        //    save = R.VY.Z; R.VY.Z = R.VZ.Z; R.VZ.Z = save;
        //  }

        //  if (diag[1] < diag[0])
        //  {
        //// swap d0 and d1
        //    save = diag[0]; diag[0] = diag[1]; diag[1] = save;

        //// swap V0 and V1
        //    save = R.VX.X; R.VX.X = R.VY.X; R.VY.X = save;
        //    save = R.VX.Y; R.VX.Y = R.VY.Y; R.VY.Y = save;
        //    save = R.VX.Z; R.VX.Z = R.VY.Z; R.VY.Z = save;
        //  }

        //  R.vy = -R.vy;
        //  R.vz = -R.vz;
        //  D.vx = vR3(diag[0],0.0,0.0);
        //  D.vy = vR3(0.0,diag[1],0.0);
        //  D.vz = vR3(0.0,0.0,diag[2]);
        //}

        #endregion

    }
}

    

